1 Experimental Design cont.

When there are further limitations as to how the experimental units can be randomly assigned to different treatment groups there are more complicated experimental designs that can be used to properly complete the analysis to reflect this limitation in randomization. We will cover Split-Plot designs and Repeated Measures.

1.1 Split-plot Design

There are situations where the logistical constraints of setting up an experiment require that there be separate random assignments of levels of factors where levels of some factor are assigned to larger experimental units (whole plots) and the levels of other factors are assigned to smaller experimental units (sub plots) within each whole plot. This unequal scale of randomization requires specific code to properly represent this experimental design in the model. While we will be using linear mixed effects models, we will wait for the mixed effects models lesson to go into the details this analysis.

1.1.1 Review the Experimental Design

The data for this example is a slightly modified version of the yield (kg/ha) trial from the RCBD example now laid out as a split-plot design (Gomez & Gomez 1984). The trial had 4 genotypes (G), 6 nitrogen levels (N or n_amount) with 3 complete replicates (rep). Within each replicate, there are six main plots (columns), and the six levels of nitrogen treatments are randomly assigned to these main plots. Within each main plot, the four genotypes are randomly assigned to subplots. Here, replicates are the ‘blocks’ in the classical split-plot terminology.

1.1.2 Importing the Data

See the “tidy” format

dat <- read_csv(here("data", "Split_Plot_Design.csv"))

## Rows: 72 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): rep, mainplot, G, N
## dbl (4): yield, row, col, n_amount
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(dat)

yield	row	col	rep	mainplot	G	N	n_amount
4520	4	1	rep1	mp1	A	N1	0
5598	2	2	rep1	mp2	A	N2	60
6192	1	3	rep1	mp3	A	N4	120
8542	2	4	rep1	mp4	A	N6	180
5806	2	5	rep1	mp5	A	N3	90
7470	1	6	rep1	mp6	A	N5	150

1.1.3 Data Description

The variables rep, mainplot, N, and G needs to be encoded as a factor instead of the character variable R uses as a default. The mutate_at function changes the selected variable (rep, mainplot, N, or G) to the chosen variable type (as.factor).

dat <- dat %>%
  mutate(across(c(rep, mainplot, N, G), as.factor)) # can use "rep:N" to include all columns between rep and N

1.1.4 Exploring the Data

We can produce the field layout of the trial using the desplot() function. The data object needs to have columns that identify the row and column of each plot in the trial. The form argument fills the color of each experimental unit by rep in this instance with a header for each rep. The text argument shows the genotype names in each plot. The remaining arguments are formatting including the main title.

desplot(data = dat,
        form = rep ~ col + row | rep, # fill color per rep, headers per rep
        text = G, cex = 1, shorten = "no", # show genotype names per plot
        col = N, # color of genotype names for each N-level
        out1 = mainplot, out1.gpar = list(col = "black"), # lines between mainplots
        out2 = row,out2.gpar = list(col = "darkgrey"), # lines between rows
                main = "Field Layout", show.key = T,key.cex = 0.7) # formatting

The field was first divided into three blocks (replicates). Then, each block was divided into six main plots (whole plot). For every block (replicates) separately, the six fertilizer treatments were randomly allocated to the main plots. Every main plot was then split into four sub plots to accommodate the four varieties. Separately for each main plot, the varieties were randomly allocated to the four sub plots. It is important to recognize that nitrogen (main plot) was randomized according to a randomized complete block design with the replicates as block. Varieties (sub plot) were also randomized according to a randomized complete block design where the main plot is the block. This type of split-plot design is the most common, but there are many other forms of the split-plot design based on how the main plot and sub plot factors are randomized

We then look at the arithmetic means and standard deviations for each genotype and nitrogen levels separately, but also their combinations. You can also arrange the tibble based on decreasing mean values

dat %>%
  group_by(G) %>%
  dplyr::summarize(mean = mean(yield),
            std.dev = sd(yield)) %>%
  arrange(desc(mean))

G	mean	std.dev
A	6553.556	1475.277
B	6155.500	1077.678
C	5563.444	1269.312
D	3642.111	1433.630

dat %>%
  group_by(N)%>%
  dplyr::summarize(mean=mean(yield),
            std.dev=sd(yield))%>%
  arrange(desc(mean))

N	mean	std.dev
N3	5866.250	832.4376
N4	5864.417	1433.5578
N5	5812.000	2349.4996
N6	5796.750	2659.6879
N2	5478.167	657.3105
N1	4054.333	671.5853

dat %>%
  group_by(N, G) %>%
  dplyr::summarize(mean = mean(yield),
            std.dev = sd(yield)) %>%
  arrange(desc(mean)) %>%
  print(n = Inf) # show more than default 10 rows

## `summarise()` has grouped output by 'N'. You can override using the `.groups`
## argument.

## # A tibble: 24 × 4
## # Groups:   N [6]
##    N     G      mean std.dev
##    <fct> <fct> <dbl>   <dbl>
##  1 N6    A     8701.   270. 
##  2 N5    A     7563.    86.9
##  3 N5    B     6951.   808. 
##  4 N4    B     6895    166. 
##  5 N4    A     6733.   490. 
##  6 N5    C     6687.   496. 
##  7 N6    B     6540.   936. 
##  8 N3    A     6400    523. 
##  9 N3    B     6259    499. 
## 10 N6    C     6065.  1097. 
## 11 N4    C     6014    515. 
## 12 N3    C     5994    101. 
## 13 N2    B     5982    684. 
## 14 N2    A     5672    458. 
## 15 N2    C     5443.   589. 
## 16 N2    D     4816    506. 
## 17 N3    D     4812    963. 
## 18 N1    D     4481.   463. 
## 19 N1    B     4306    646. 
## 20 N1    A     4253.   248. 
## 21 N4    D     3816   1311. 
## 22 N1    C     3177.   453. 
## 23 N5    D     2047.   703. 
## 24 N6    D     1881.   407.

A simple plot of the raw data also helps to visualize the distribution of the yield for each genotype nitrogen treatment combination. The point geometry produces a scatter plot. The ylim argument forced the y-axis to start at 0. The theme_classic argument produces a clearer plot format.

ggplot(data = dat,
       aes(y = yield, x = N, color = N)) + # different colors for the level of nitrogen
  facet_wrap(~ G) + #facet per G level - "labeller = label_both" adds the variable name to each facet graph
  geom_point() +
  scale_y_continuous(limits = c(0, NA), # make y = axis start at 0
                     expand = expansion(mult = c(0, 0.1)) # no space below 0
                     ) +
  scale_x_discrete(name = NULL) + # x-axis with no name label
  theme_bw() + # alternative plot theme with more grid lines - personal preference for theme_classic
  theme(legend.position = "bottom") # legend on bottom

1.1.5 Fitting ANOVA Models with R

We fit a linear model with yield as the response variable. The remaining parts of the model can be broken up into two parts - the design effects and treatment effects. The treatment effects are genotype, nitrogen level, and their interaction as the fixed effects. For the split-plot design there are two randomization units to be represented in the linear model - the main plots and sub plots. Each randomization unit needs to be represented by a random effect so each randomization unit has its own error term.

The goal is to determine if there is a difference in yield between varieties and where those differences are.

The interaction notation N * G is short for N + G + N:G

(1|rep) represents the block level variance - random effect (1|rep:mainplot) represents the whole-plot (main plot) error - random effect The residual represents the sub-plot variability not captured by mainplot variance

mod_re <- lmerTest::lmer(yield ~ N * G + (1|rep) + (1|rep:mainplot), data = dat)

## boundary (singular) fit: see help('isSingular')

summary(mod_re) # We can ignore the Singular fit warning message for now

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: yield ~ N * G + (1 | rep) + (1 | rep:mainplot)
##    Data: dat
## 
## REML criterion at convergence: 780.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8506 -0.5356 -0.1311  0.4778  1.9272 
## 
## Random effects:
##  Groups       Name        Variance     Std.Dev. 
##  rep:mainplot (Intercept)  50855.53547 225.51172
##  rep          (Intercept)      0.00307   0.05541
##  Residual                 349441.51914 591.13579
## Number of obs: 72, groups:  rep:mainplot, 18; rep, 3
## 
## Fixed effects:
##             Estimate Std. Error       df t value            Pr(>|t|)    
## (Intercept)  4252.67     365.28    45.78  11.642 0.00000000000000279 ***
## NN2          1419.33     516.59    45.78   2.748             0.00856 ** 
## NN3          2147.33     516.59    45.78   4.157             0.00014 ***
## NN4          2480.00     516.59    45.78   4.801 0.00001725039156940 ***
## NN5          3310.67     516.59    45.78   6.409 0.00000007184143597 ***
## NN6          4448.00     516.59    45.78   8.610 0.00000000003937369 ***
## GB             53.33     482.66    36.00   0.110             0.91263    
## GC          -1075.33     482.66    36.00  -2.228             0.03222 *  
## GD            228.67     482.66    36.00   0.474             0.63853    
## NN2:GB        256.67     682.58    36.00   0.376             0.70911    
## NN3:GB       -194.33     682.58    36.00  -0.285             0.77750    
## NN4:GB        109.00     682.58    36.00   0.160             0.87402    
## NN5:GB       -666.00     682.58    36.00  -0.976             0.33572    
## NN6:GB      -2213.67     682.58    36.00  -3.243             0.00255 ** 
## NN2:GC        846.00     682.58    36.00   1.239             0.22321    
## NN3:GC        669.33     682.58    36.00   0.981             0.33334    
## NN4:GC        356.67     682.58    36.00   0.523             0.60451    
## NN5:GC        199.33     682.58    36.00   0.292             0.77194    
## NN6:GC      -1560.00     682.58    36.00  -2.285             0.02828 *  
## NN2:GD      -1084.67     682.58    36.00  -1.589             0.12079    
## NN3:GD      -1816.67     682.58    36.00  -2.661             0.01155 *  
## NN4:GD      -3145.33     682.58    36.00  -4.608 0.00004945084378565 ***
## NN5:GD      -5745.33     682.58    36.00  -8.417 0.00000000050133641 ***
## NN6:GD      -7048.67     682.58    36.00 -10.326 0.00000000000261524 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

mod_re %>% car::Anova(type = 3, test.statistic = "F")

	F	Df	Df.res	Pr(>F)
(Intercept)	135.538147	1	45.78314	0.0000000
N	17.621854	5	41.96547	0.0000000
G	3.016614	3	36.00000	0.0424007
N:G	13.235986	15	36.00000	0.0000000

epsilon <- residuals(mod_re)

1.1.6 Sums of Squares Review

There are three types of sums of squares (1, 2, or 3). This is an important consideration when the data is unbalanced (missing data). To demonstrate the differences we will consider a model that includes two categorical independent variables (Factor A and Factor B). The model will test the two main effects and the interaction between the two factors. Type 1 sums of squares is also called “sequential” sums of square that tests for Factor A, then Factor B, then the interaction of Factor A and B. When the data is unbalanced the sums of squares will give a different result depending on which main effect is considered first. It is rare that you want to test the main effects sequentially, so normally you should use Type 2 or 3 sums of squares.

Type 2 sums of squares is used only when you are not testing an interaction, and tests each main effect after the other main effect. Type 3 sums of squares tests each main effect after accounting for all other effects, including interactions;it is commonly used when you are testing an interaction. If you test the interaction and it is not significant, then you can test the main effects alone and use the type 2 sums of squares - it is more powerful than type 3.

This link can provide more details for your reference.

1.1.7 Assumptions

Multiple ways to assess model assumptions - the autoplot function used before doesn’t work with lmer models.

qqnorm(epsilon)
qqline(epsilon)

plot(mod_re, type = c("p", "smooth"), col.line = 1)

plot(residuals(mod_re) ~ fitted(mod_re))
abline(h = 0)

lattice::qqmath(mod_re)

check_model(mod_re, check = c("normality", "homogeneity"))

shapiro.test(epsilon)

## 
##  Shapiro-Wilk normality test
## 
## data:  epsilon
## W = 0.96729, p-value = 0.05772

dat$resid <- residuals(mod_re)
car::leveneTest(resid ~ N * G, data = dat)

	Df	F value	Pr(>F)
group	23	0.4656747	0.9752327
	48	NA	NA

Results indicate the satisfaction of linear model assumptions, so we can move onto multiple comparisons tests.

1.1.8 Mean Comparisons

Following a significant F-test, we compare the variety means. We will show all pairwise comparisons using the colon between the two factors.

all_mean_comparisons_tukey_re <- mod_re %>%
  emmeans(specs = ~ N:G) %>%
  cld(Letters = letters)
all_mean_comparisons_tukey_re

	N	G	emmean	SE	df	lower.CL	upper.CL	.group
24	N6	D	1880.667	365.2839	45.78314	1145.294	2616.039	a
23	N5	D	2046.667	365.2839	45.78314	1311.294	2782.039	a
13	N1	C	3177.333	365.2839	45.78314	2441.961	3912.706	ab
22	N4	D	3816.000	365.2839	45.78314	3080.628	4551.372	abc
1	N1	A	4252.667	365.2839	45.78314	3517.294	4988.039	bcd
7	N1	B	4306.000	365.2839	45.78314	3570.628	5041.372	bcd
19	N1	D	4481.333	365.2839	45.78314	3745.961	5216.706	bcde
21	N3	D	4812.000	365.2839	45.78314	4076.628	5547.372	bcdef
20	N2	D	4816.000	365.2839	45.78314	4080.628	5551.372	bcdef
14	N2	C	5442.667	365.2839	45.78314	4707.294	6178.039	cdefg
2	N2	A	5672.000	365.2839	45.78314	4936.628	6407.372	cdefgh
8	N2	B	5982.000	365.2839	45.78314	5246.628	6717.372	defgh
15	N3	C	5994.000	365.2839	45.78314	5258.628	6729.372	defgh
16	N4	C	6014.000	365.2839	45.78314	5278.628	6749.372	defgh
18	N6	C	6065.333	365.2839	45.78314	5329.961	6800.706	defgh
9	N3	B	6259.000	365.2839	45.78314	5523.628	6994.372	defgh
3	N3	A	6400.000	365.2839	45.78314	5664.628	7135.372	efgh
12	N6	B	6540.333	365.2839	45.78314	5804.961	7275.706	fgh
17	N5	C	6687.333	365.2839	45.78314	5951.961	7422.706	fgh
4	N4	A	6732.667	365.2839	45.78314	5997.294	7468.039	fghi
10	N4	B	6895.000	365.2839	45.78314	6159.628	7630.372	ghi
11	N5	B	6950.667	365.2839	45.78314	6215.294	7686.039	ghi
5	N5	A	7563.333	365.2839	45.78314	6827.961	8298.706	hi
6	N6	A	8700.667	365.2839	45.78314	7965.294	9436.039	i

1.1.9 Data Visualization

Create a plot that displays both the raw data and results from the multiple comparisons tests of the adjusted means from the linear model. First create a tibble with an additional column that combined the N and G variables, and then sorts the rows on increasing values of emmean. Then create a similar table for the raw data.

all_mean_comparisons_tukey_re <- all_mean_comparisons_tukey_re %>%
  as_tibble() %>%
  mutate(N_G = paste0(N,"-",G)) %>%
  mutate(N_G = fct_reorder(N_G, emmean))

dat <- dat %>%
  mutate(N_G = paste0(N,"-",G)) %>%
  mutate(N_G = fct_relevel(N_G, levels(all_mean_comparisons_tukey_re$N_G)))

RCBD_Plot_tukey <- ggplot() +
  geom_point(data = dat, aes(y = yield, x = N_G, color = N), position = position_nudge(x = -0.4)) + # black dots representing the raw data
  geom_boxplot(data =dat, aes(y = yield, x = N_G), width = 0.25, outlier.shape = NA, position = position_nudge(x = -0.2)) + # black boxplot
  geom_point(data = all_mean_comparisons_tukey_re, aes(y = emmean, x = N_G), size = 2, color = "red") + #  red dots representing the adjusted means
  geom_errorbar(data = all_mean_comparisons_tukey_re, aes(ymin = lower.CL, ymax = upper.CL, x = N_G), color = "red", width = 0.1) + # red error bars representing the confidence limits of the adjusted means
  geom_text(data = all_mean_comparisons_tukey_re, aes(y = 10000, x = N_G, label = str_trim(.group)), color = "red", position = position_nudge(x = -0.2), hjust = 0, angle = 90) + # red letters 
  scale_y_continuous(name = "Yield", limits = c(0, NA), expand = expansion(mult = c(0, 0.1))) +
  scale_x_discrete(name = "Nitrogen-Genotype combination") +
  theme_classic() + # clearer plot format 
labs(caption=str_wrap("Black dots represent raw data. Red dots and error bars represent estimated marginal means ± 95% confidence interval per group. Means not sharing any letter are significantly different by the t-test at the 5% level of significance.", width = 120)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), legend.position = "bottom")
RCBD_Plot_tukey

1.2 Repeated Measures

When repeated observations are made on the same individual experimental unit this violates one of the linear model assumptions - no autocorrelation among errors. This analysis requires some specific steps to account for the temporal autocorrelation. We will first fit a model that assumes no autocorrelation. Then we will fit multiple models that assume different types of temporal autocorrelation and compare them to see which fits the data best.

1.2.1 Review the Experimental Design

Data from a sorghum trial laid out as a randomized complete block design with five blocks and four sorghum varieties as the only treatment factor. The response variable - leaf area index - was measured in five consecutive weeks on each plot starting two weeks after emergence. The repeated observation variable (week) needs to be treated as a factor - all measurements were made at discrete time points. This dataset is comprised of 100 values resembling longitudinal data - repeated measures or time series analysis.

The week factor is not a treatment factor that can be randomized. The repeated measurements of the same experimental unit are likely to be serially correlated, which needs to be accounted for to produce a non-biased estimate of the differences between varieties.

1.2.2 Importing the Data

We also include the steps to format the necessary variables to be factors.

# data - reformatting agriTutorial::sorghum
dat <- agriTutorial::sorghum %>% # data from agriTutorial package
  rename(block = Replicate, 
         weekF = factweek,  # week as factor
         weekN = varweek,   # week as numeric/integer
         plot  = factplot) %>% 
  mutate(variety = paste0("var", variety),    # variety id
         block   = paste0("block", block),    # block id
         weekF   = paste0("week", weekF),     # week id
         plot    = paste0("plot", plot),      # plot id
         unit    = paste0("obs", 1:n() )) %>% # observation id - needed for glmmTMB model
  mutate_at(vars(variety:plot, unit), as.factor)
head(dat)

y	variety	block	weekF	plot	weekN	varblock	unit
5.00	var1	block1	week1	plot1	1	1	obs1
4.84	var1	block1	week2	plot1	2	1	obs2
4.02	var1	block1	week3	plot1	3	1	obs3
3.75	var1	block1	week4	plot1	4	1	obs4
3.13	var1	block1	week5	plot1	5	1	obs5
4.42	var1	block2	week1	plot2	1	2	obs6

1.2.3 Exploring the Data

The lack of row and column variables prevent us from using the desplot function to create a figure that shows the experimental design. Descriptive tables will show the arithmetic means and standard deviations for leaf area index per variety. We can use the pivot_wider function to get mean values for each week.

dat %>%
  group_by(variety) %>% 
  summarize(mean = mean(y, na.rm = TRUE), std.dev = sd(y, na.rm = TRUE)) %>% 
  arrange(desc(mean)) # sort

variety	mean	std.dev
var3	4.6340	0.7568741
var4	4.6052	0.8024033
var2	4.5900	0.7634352
var1	3.4960	0.7601425

dat %>%
  group_by(weekF, variety) %>% 
  summarize(mean = mean(y, na.rm = TRUE)) %>% 
  pivot_wider(names_from = weekF, values_from = mean) # pivot wider creates the table with each week as a column

## `summarise()` has grouped output by 'weekF'. You can override using the
## `.groups` argument.

variety	week1	week2	week3	week4	week5
var1	4.242	4.046	3.406	3.086	2.700
var2	5.148	4.980	4.594	4.220	4.008
var3	4.888	5.074	4.714	4.476	4.018
var4	5.154	4.944	4.706	4.388	3.834

We can create an animated plot that shows the change in leaf area index by variety across the five weeks.

var_colors <- c("#8cb369", "#f4a259", "#5b8e7d", "#bc4b51")
names(var_colors) <- dat$variety %>% levels()

gganimate_plot <- ggplot(
  data = dat, aes(y = y, x = weekF,
                  group = variety,
                  color = variety)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(alpha = 0.5, size = 3) +
  scale_y_continuous(
    name = "Leaf area index",
    limits = c(0, 6.5),
    expand = c(0, 0),
    breaks = c(0:6)
  ) +
  scale_color_manual(values = var_colors) +
  theme_bw() +
  theme(legend.position = "bottom", 
        axis.title.x = element_blank()) +
  transition_time(weekN) +
  shadow_mark(exclude_layer = 2) 

animate(gganimate_plot, renderer = gifski_renderer()) # render gif

1.2.4 Fitting ANOVA Models with R

An inefficient approach is to model a subset of the data from one week using either a fixed effect for block or random effect.

dat.wk1 <- dat %>% filter(weekF == "week1") # subset data from first week only
mod.wk1 <- lm(y ~ variety + block, data = dat.wk1)
anova(mod.wk1)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
variety	3	2.76036	0.9201200	46.59032	0.0000007
block	4	8.31897	2.0797425	105.30786	0.0000000
Residuals	12	0.23699	0.0197492	NA	NA

mod.wk1_re <- lmerTest::lmer(y ~ variety + (1|block), data = dat.wk1)
anova(mod.wk1_re)

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
variety	2.76036	0.92012	3	12	46.59032	0.0000007

dat.wk2 <- dat %>% filter(weekF == "week2")
mod.wk2 <- lm(y ~ variety + block, data = dat.wk2)
anova(mod.wk2)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
variety	3	3.45322	1.151073	55.97245	0.0000003
block	4	8.60998	2.152495	104.66788	0.0000000
Residuals	12	0.24678	0.020565	NA	NA

mod.wk2_re <- lmerTest::lmer(y ~ variety + (1|block), data = dat.wk2)
anova(mod.wk2_re)

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
variety	3.45322	1.151073	3	12	55.97245	0.0000003

#Repeat with other weeks

We can also fit a model with all weeks combined.

It is reasonable to assume that the treatment effects evolve over time and thus are week-specific. Furthermore, there could be variability among blocks over the weeks. The glmmTMB function allows us more flexibility than the lmer function.

mod.iid <- glmmTMB(y ~ weekF * variety + (1|block) + (1|unit), # add random unit term to mimic error variance
                 dispformula = ~ 0, # fix original error variance to 0
                 REML = TRUE, # needs to be stated since default is ML
                 data = dat)

mod.iid %>% broom.mixed::tidy(effects = "ran_pars", scales = "vcov")

effect	component	group	term	estimate
ran_pars	cond	block	var__(Intercept)	0.4144806
ran_pars	cond	unit	var__(Intercept)	0.0271694

These results assume independent and homogeneous error term, which MAY be wrong. We will now test for possible covariance structures for the error term that improve overall model fit. This should be done before assessing the significance of any effect - the model assumptions should be met before assessing the significant differences among the varieties.

1.2.5 Autocorrelated Errors

Experimental units on which repeated observations are made are often referred to as subjects. We would now like to allow measurements taken on the same subjects (plot in this case) to be serially correlated, while observations on different subjects are still considered independent. More specifically, we want the errors of the respective observations to be correlated in our model. It is also reasonable to assume a weaker correlation between errors that are further apart in time between measurements.

# Function to create correlation matrices
make_corr_matrix <- function(type, n=5, rho=0.6) {
  mat <- matrix(0, n, n)
  if(type=="AR1") {
    for(i in 1:n) for(j in 1:n) mat[i,j] <- rho^abs(i-j)
  } else if(type=="CS") {
    mat <- matrix(rho, n, n)
    diag(mat) <- 1
  } else if(type=="Toeplitz") {
    rho_seq <- seq(rho, rho^n, length.out=n)
    for(i in 1:n) for(j in 1:n) mat[i,j] <- rho_seq[abs(i-j)+1]
    diag(mat) <- 1
  } else if(type=="UN") {
    mat <- matrix(NA, n, n) # unstructured: all correlations potentially different
  }
  rownames(mat) <- colnames(mat) <- paste0("t", 1:n)
  mat
}

plot_corr_matrix <- function(mat, title) {
  melted <- melt(mat, na.rm=TRUE)
  ggplot(melted, aes(Var1, Var2, fill=value)) +
    geom_tile(color="white") +
    scale_fill_gradient2(low="blue", high="red", mid="white", midpoint=0.5, limits=c(0,1)) +
    theme_minimal() + coord_fixed() +
    ggtitle(title) + xlab("") + ylab("") +
    theme(axis.text.x = element_text(angle=45, hjust=1))
}

# Plot examples
plot_corr_matrix(make_corr_matrix("AR1"), "AR(1) Correlation")

plot_corr_matrix(make_corr_matrix("CS"), "Compound Symmetry (CS) Correlation")

plot_corr_matrix(make_corr_matrix("Toeplitz"), "Toeplitz Correlation")

plot_corr_matrix(make_corr_matrix("UN"), "Unstructured (UN) Correlation")

AR(1) decreases with lag
CS is constant
Toeplitz varies arbitrarily by lag
UN allows every entry to differ

1.2.6 AR(1)

The most popular correlation structure is called first order autoregressive AR(1). This error structure works best when all time points are equally spaced. Other variance structures are possible and should also be considered.

mod.AR1 <- glmmTMB(formula = y ~ weekF * variety + (1|block) +
                   ar1(weekF + 0 | plot), # ar1 structure as random term to mimic error var, and + 0 is used to remove intercept term from random term covariance specification, which is required to properly define the covariance matrix for repeated measures
                   dispformula = ~ 0, # fix original error variance to 0
                   REML = TRUE,       # needs to be stated since default = ML
                   data = dat) 

# Extract variance component estimates
mod.AR1 %>% broom.mixed::tidy(effects = "ran_pars", scales = "vcov")

effect	component	group	term	estimate
ran_pars	cond	block	var__(Intercept)	0.3916641
ran_pars	cond	plot	var__weekFweek1	0.0325258
ran_pars	cond	plot	var__weekFweek2	0.0325258
ran_pars	cond	plot	var__weekFweek3	0.0325258
ran_pars	cond	plot	var__weekFweek4	0.0325258
ran_pars	cond	plot	var__weekFweek5	0.0325258
ran_pars	cond	plot	cov__weekFweek1.weekFweek2	0.0243886
ran_pars	cond	plot	cov__weekFweek1.weekFweek3	0.0182872
ran_pars	cond	plot	cov__weekFweek1.weekFweek4	0.0137122
ran_pars	cond	plot	cov__weekFweek1.weekFweek5	0.0102817
ran_pars	cond	plot	cov__weekFweek2.weekFweek3	0.0243886
ran_pars	cond	plot	cov__weekFweek2.weekFweek4	0.0182872
ran_pars	cond	plot	cov__weekFweek2.weekFweek5	0.0137122
ran_pars	cond	plot	cov__weekFweek3.weekFweek4	0.0243886
ran_pars	cond	plot	cov__weekFweek3.weekFweek5	0.0182872
ran_pars	cond	plot	cov__weekFweek4.weekFweek5	0.0243886

1.2.7 Compound Symmetry

Compound symmetry can be modeled as standard - homogeneous (nlme only - not showing here) or heterogeneous (glmmTMB) to model the error term. Assumes the correlation is constant regardless of how far apart the measurements are. The heterogeneous compound symmetry is a simple extension of the homogeneous compound symmetry structure with more parameters to allow for the variances along the diagonal of the matrix to be different.

mod.hCS <- glmmTMB(formula = y ~ weekF * variety + (1|block) +
                   cs(weekF + 0 | plot), # hcs structure as random term to mimic error var
                   dispformula = ~ 0, # fix original error variance to 0
                   REML = TRUE,       # needs to be stated since default = ML
                   data = dat) 

# show variance components
mod.hCS %>% broom.mixed::tidy(effects = "ran_pars", scales = "vcov")

effect	component	group	term	estimate
ran_pars	cond	block	var__(Intercept)	0.8033429
ran_pars	cond	plot	var__weekFweek1	0.0552712
ran_pars	cond	plot	var__weekFweek2	0.0584155
ran_pars	cond	plot	var__weekFweek3	0.1155951
ran_pars	cond	plot	var__weekFweek4	0.1247201
ran_pars	cond	plot	var__weekFweek5	0.1990451
ran_pars	cond	plot	cov__weekFweek1.weekFweek2	0.0522395
ran_pars	cond	plot	cov__weekFweek1.weekFweek3	0.0734860
ran_pars	cond	plot	cov__weekFweek1.weekFweek4	0.0763314
ran_pars	cond	plot	cov__weekFweek1.weekFweek5	0.0964297
ran_pars	cond	plot	cov__weekFweek2.weekFweek3	0.0755474
ran_pars	cond	plot	cov__weekFweek2.weekFweek4	0.0784726
ran_pars	cond	plot	cov__weekFweek2.weekFweek5	0.0991347
ran_pars	cond	plot	cov__weekFweek3.weekFweek4	0.1103885
ran_pars	cond	plot	cov__weekFweek3.weekFweek5	0.1394541
ran_pars	cond	plot	cov__weekFweek4.weekFweek5	0.1448538

1.2.8 Toeplitz

Toeplitz (glmmTMB only) can be calculated faster than similar AR(1) approach. Toeplitz is unique from AR(1) in that the correlations do not necessarily have the same pattern as is the case in AR(1). This correlation structure is only possible using the glmmTMB function.

mod.Toep <- glmmTMB(formula = y ~ weekF * variety + (1|block) +
                   toep(weekF + 0 | plot), # toep structure as random term to mimic err var
                   dispformula = ~ 0, # fix original error variance to 0
                   REML = TRUE,       # needs to be stated since default = ML
                   data = dat)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

# show variance components
mod.Toep %>% broom.mixed::tidy(effects = "ran_pars", scales = "vcov")

effect	component	group	term	estimate
ran_pars	cond	block	var__(Intercept)	0.3436258
ran_pars	cond	plot	var__weekFweek1	0.0393570
ran_pars	cond	plot	var__weekFweek2	0.0392134
ran_pars	cond	plot	var__weekFweek3	0.0236161
ran_pars	cond	plot	var__weekFweek4	0.0342313
ran_pars	cond	plot	var__weekFweek5	0.0301218
ran_pars	cond	plot	cov__weekFweek1.weekFweek2	0.0300907
ran_pars	cond	plot	cov__weekFweek1.weekFweek3	0.0175201
ran_pars	cond	plot	cov__weekFweek1.weekFweek4	0.0099941
ran_pars	cond	plot	cov__weekFweek1.weekFweek5	-0.0005475
ran_pars	cond	plot	cov__weekFweek2.weekFweek3	0.0233091
ran_pars	cond	plot	cov__weekFweek2.weekFweek4	0.0210548
ran_pars	cond	plot	cov__weekFweek2.weekFweek5	0.0093579
ran_pars	cond	plot	cov__weekFweek3.weekFweek4	0.0217781
ran_pars	cond	plot	cov__weekFweek3.weekFweek5	0.0153273
ran_pars	cond	plot	cov__weekFweek4.weekFweek5	0.0245955

1.2.9 Unstructured Variance

The unstructured variance structure is the most “liberal” that allows every term to be different. In doing so it requires fitting the most parameters.

mod.UN <- glmmTMB(formula = y ~ weekF * variety + (1|block) +
                   us(weekF + 0 | plot), # us structure as random term to mimic error var
                   dispformula = ~ 0, # fix original error variance to 0
                   REML = TRUE,       # needs to be stated since default = ML
                   data = dat) 

# show variance components
mod.UN %>% broom.mixed::tidy(effects = "ran_pars", scales = "vcov")

effect	component	group	term	estimate
ran_pars	cond	block	var__(Intercept)	0.4181621
ran_pars	cond	plot	var__weekFweek1	0.0243949
ran_pars	cond	plot	var__weekFweek2	0.0264585
ran_pars	cond	plot	var__weekFweek3	0.0223323
ran_pars	cond	plot	var__weekFweek4	0.0312693
ran_pars	cond	plot	var__weekFweek5	0.0408654
ran_pars	cond	plot	cov__weekFweek1.weekFweek2	0.0199855
ran_pars	cond	plot	cov__weekFweek1.weekFweek3	0.0087999
ran_pars	cond	plot	cov__weekFweek1.weekFweek4	0.0126196
ran_pars	cond	plot	cov__weekFweek1.weekFweek5	-0.0028035
ran_pars	cond	plot	cov__weekFweek2.weekFweek3	0.0109330
ran_pars	cond	plot	cov__weekFweek2.weekFweek4	0.0153952
ran_pars	cond	plot	cov__weekFweek2.weekFweek5	-0.0009543
ran_pars	cond	plot	cov__weekFweek3.weekFweek4	0.0248921
ran_pars	cond	plot	cov__weekFweek3.weekFweek5	0.0197851
ran_pars	cond	plot	cov__weekFweek4.weekFweek5	0.0246111

1.2.10 Model Selection

In order to select the best model here, we can simply compare their AIC values, since all models are identical regarding their fixed effects. The smaller the value of AIC, the better the fit.

AICcmodavg::aictab(
  cand.set = list(mod.iid, mod.hCS, mod.AR1, mod.Toep, mod.UN), 
  modnames = c("iid", "hCS", "AR1", "Toeplitz", "UN"),
  second.ord = FALSE,   # get AIC instead of AICc
  sort = TRUE
)

	Modnames	K	AIC	Delta_AIC	ModelLik	AICWt	LL	Cum.Wt
5	UN	36	-9.7867784	0.000000	1.0000000	0.9335234	40.893389	0.9335234
3	AR1	23	-4.3473941	5.439384	0.0658950	0.0615146	25.173697	0.9950379
4	Toeplitz	30	0.6958276	10.482606	0.0052934	0.0049415	29.652086	0.9999794
2	hCS	27	11.6540463	21.440825	0.0000221	0.0000206	21.172977	1.0000000
1	iid	22	37.6615506	47.448329	0.0000000	0.0000000	3.169225	1.0000000

According to AIC values the UN correlation structure has the lowest AIC value as no others are within 2. It is common to use Delta AIC < 2 to select a model. We will now review the UN model in the following steps.

The Anova function produces the table of F-Statistics to determine if there significant differences among the different levels of the varieties tested.

mod.UN %>% car::Anova(type = 3) #type II sums of squares is not appropriate with interaction

	Chisq	Df
(Intercept)	203.3025	1
weekF	299.8838	4
variety	113.1530	3
weekF:variety	111.1423	12

epsilon <- residuals(mod.UN) #used to test model assumptions

1.2.11 Assumptions

Before we look at any multiple comparisons we need to confirm that the model assumptions have been met.

qqnorm(epsilon)
qqline(epsilon)

plot(residuals(mod.UN) ~ fitted(mod.UN))
abline(h = 0)

#compared to iid model
plot(residuals(mod.iid) ~ fitted(mod.iid)) #less obvious 
abline(h = 0)

check_model(mod.UN, check = c("normality", "homogeneity"))

check_model(mod.iid, check = c("normality", "homogeneity"))

shapiro.test(epsilon)

## 
##  Shapiro-Wilk normality test
## 
## data:  epsilon
## W = 0.9935, p-value = 0.9162

car::leveneTest(epsilon ~ weekF * variety, data = dat)

	Df	F value	Pr(>F)
group	19	1.595601	0.0776386
	80	NA	NA

Results indicate the satisfaction of linear model assumptions, so we can move onto multiple comparisons tests.

1.2.12 Mean Comparisons

all_mean_comparisons_tukey_UN <- mod.UN %>% # all possible comparisons aren't needed based on primary objectives for comparing varieties within week
  emmeans(specs = ~ variety:weekF) %>%
  cld(Letters = letters)

withinWeek_mean_comparisons_tukey_UN <- mod.UN %>% 
  emmeans(specs = ~ variety|weekF) %>%
  cld(Letters = letters)

1.2.13 Data Visualization

Create a plot that displays both the raw data and results from the multiple comparisons tests of the adjusted means from the linear model. First create a tibble with an additional column that combined the week and variety variables, and then sorts the rows on increasing values of emmean. Then create a similar table for the raw data.

withinWeek_UN_Plot_tukey <- ggplot() +
  facet_wrap(~ weekF, labeller = label_both, nrow = 1) + # facet per G level
  geom_point(data = dat, aes(y = y, x = variety)) + # dots representing the raw data
  geom_point(data = withinWeek_mean_comparisons_tukey_UN, aes(y = emmean, x = variety), color = "red", position = position_nudge(x = 0.1)) + # red dots representing the adjusted means
  geom_errorbar(data = withinWeek_mean_comparisons_tukey_UN, aes(ymin = lower.CL, ymax = upper.CL, x = variety), color = "red", width = 0.1, position = position_nudge(x = 0.1)) + # red error bars representing the confidence limits of the adjusted means
  geom_text(data = withinWeek_mean_comparisons_tukey_UN, aes(y = emmean, x = variety, label = str_trim(.group)), color = "red", position=position_nudge(x = 0.2), hjust = 0) + # red letters 
  scale_y_continuous(name = "Yield", limits = c(0, NA), expand = expansion(mult = c(0, 0.1))) +
  scale_x_discrete(name = NULL) +
  theme_classic() + # clearer plot format 
labs(caption = str_wrap("Black dots represent raw data. Red dots and error bars represent estimated marginal means ± 95% confidence interval per group. Means not sharing any letter are significantly different by the t-test at the 5% level of significance.", width = 120)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), legend.position = "bottom")
withinWeek_UN_Plot_tukey

1.3 Practice Example

We will use the data from an experiment conducted in a greenhouse in Silver Bay, Minnesota. Plants exposed to different (three) treatments were observed over a 40 day period. Make sure to modify the variable types to factor for the repeated measures analysis. Fit a linear model assuming no random effects, a linear mixed effects model with a random intercept for each plant, and a linear mixed effects model that tests unstructured variance, heterogeneous compound symmetry, autoregressive order-1, homogeneous diagonal variance, and heterogeneous diagonal variance structures.

library(agridat)
data(pederson.lettuce.repeated)
dat <- pederson.lettuce.repeated
library(lattice)
dat <- dat[order(dat$day),]
xyplot(weight ~ day|trt, dat, type = 'l', group = plant, layout = c(3, 1),
main = "pederson.lettuce.repeated")

1.3.1 Tips if you get stuck

You really want to see the answers without trying yourself???

# Remember to change the necessary variables to factors
str(dat)

## 'data.frame':    594 obs. of  4 variables:
##  $ plant : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ day   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ trt   : Factor w/ 3 levels "T1","T2","T3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ weight: num  44 48.4 55.5 48.6 48.1 42.6 55.4 47.5 48.8 51.2 ...

dat$plant <- as.factor(dat$plant)
dat$day <- as.factor(dat$day)

# Fitting linear model with no random effect
mod_lm <- lm(weight ~ trt * day, data = dat)

# Fitting linear mixed model with random effect
mod <- glmmTMB(weight ~ trt * day + (1|plant), data = dat)

# Code to try when issues of convergence pop up
# control = glmmTMBControl(optimizer = optim, optArgs = list(method = "BFGS"))
# control = glmmTMBControl(maxfun = 20000)

mod.UN <- glmmTMB(weight ~ trt * day + us(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat)

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')

mod.UN <- glmmTMB(weight ~ trt * day + us(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat,
                control = glmmTMBControl(optimizer = optim, optArgs = list(method = "BFGS")))

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; . See vignette('troubleshooting'), help('diagnose')

# Unstructured Correlation - fit estimate for each pair of observations over time
# Overparameterized model warning message - makes sense this unstructured covariance option fits a covariance estimate for every pair of observations


mod.CS <- glmmTMB(weight ~ trt * day + cs(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat)

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')

mod.CS <- glmmTMB(weight ~ trt * day + cs(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat,
                control = glmmTMBControl(optimizer = optim, optArgs = list(method = "BFGS")))
mod.HomD <- glmmTMB(weight ~ trt * day + homdiag(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat)
mod.HetD <- glmmTMB(weight ~ trt * day + diag(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat)

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')

# Convergence warning - try different optimizer
# Very common to get convergence warnings for small datasets of some covariance structures
mod.HetD <- glmmTMB(weight ~ trt * day + diag(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat,
                control = glmmTMBControl(optimizer = optim, optArgs = list(method = "BFGS")))
mod.AR1 <- glmmTMB(weight ~ trt * day + ar1(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat) 

AICcmodavg::aictab(
  cand.set = list(mod, mod.CS, mod.HomD, mod.HetD, mod.AR1), 
  modnames = c("iid", "hCS", "HomD", "HetD", "Ar1"),
  second.ord = FALSE) # get AIC instead of AICc

	Modnames	K	AIC	Delta_AIC	ModelLik	AICWt	LL	Cum.Wt
5	Ar1	35	2257.262	0.0000	1	1	-1093.631	1
2	hCS	45	2718.281	461.0183	0	0	-1314.140	1
1	iid	35	2999.319	742.0571	0	0	-1464.660	1
4	HetD	44	3409.611	1152.3487	0	0	-1660.806	1
3	HomD	34	3498.414	1241.1513	0	0	-1715.207	1

# AR1 model is best fitting model
# Now to check assumptions and make any transformations as needed
mod.AR1 <- glmmTMB(weight ~ trt * day + ar1(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat) 
plot(residuals(mod.AR1) ~ fitted(mod.AR1)) # less obvious 
abline(h = 0)

check_model(mod.AR1, check = c("normality", "homogeneity"))

qqnorm(residuals(mod.AR1))
qqline(residuals(mod.AR1))

shapiro.test(residuals(mod.AR1))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod.AR1)
## W = 0.93534, p-value = 0.000000000000002314

# Try polynomial term for day - orthogonal polynomial (not raw polynomial coefficients)
mod.AR1_poly <- glmmTMB(weight ~ trt * poly(day, 2) + ar1(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat) 
qqnorm(residuals(mod.AR1_poly))
qqline(residuals(mod.AR1_poly))

shapiro.test(residuals(mod.AR1_poly))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod.AR1_poly)
## W = 0.97337, p-value = 0.000000006418

plot(residuals(mod.AR1_poly) ~ fitted(mod.AR1_poly)) # less obvious 
abline(h = 0)

# Other resources to assess assumptions for more complicated models
sims = simulateResiduals(mod.AR1) 
plot(sims, quantreg = T)

sims = simulateResiduals(mod.AR1_poly) 
plot(sims, quantreg = T)

# KS Test: expected to follow a uniform distribution - not necessarily normal distribution values close to 1.0 indicate a good fit while values close to 0.0 indicate a poor fit.

mod.AR1 <- glmmTMB(log(weight) ~ trt * day + ar1(day + 0|plant),
                dispformula = ~ 0, REML = TRUE, data = dat) 

sims = simulateResiduals(mod.AR1) 
plot(sims, quantreg = T)

qqnorm(residuals(mod.AR1))
qqline(residuals(mod.AR1))

shapiro.test(residuals(mod.AR1))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod.AR1)
## W = 0.92093, p-value < 0.00000000000000022

plot(residuals(mod.AR1) ~ fitted(mod.AR1)) # less obvious 
abline(h = 0)

car::Anova(mod.AR1, type = "III")

	Chisq	Df	Pr(>Chisq)
(Intercept)	56476.784784	1	0.0000000
trt	1.069904	2	0.5856973
day	2123.165480	10	0.0000000
trt:day	225.952277	20	0.0000000

within_means_tukey <- mod.AR1 %>%
  emmeans(~ trt|day) %>%
  cld(Letters = letters) %>%
  as_tibble() %>%
  rename(weight = emmean)
within_means_tukey

trt	day	weight	SE	df	lower.CL	upper.CL	.group
T2	1	3.886148	0.0163905	559	3.853953	3.918342	a
T1	1	3.895169	0.0163905	559	3.862974	3.927363	a
T3	1	3.909896	0.0163905	559	3.877702	3.942091	a
T1	5	3.909989	0.0163905	559	3.877794	3.942183	a
T2	5	3.915754	0.0163905	559	3.883559	3.947948	a
T3	5	3.954929	0.0163905	559	3.922734	3.987123	a
T1	8	3.945891	0.0163905	559	3.913696	3.978085	a
T2	8	3.947583	0.0163905	559	3.915388	3.979777	a
T3	8	3.974883	0.0163905	559	3.942689	4.007078	a
T2	12	3.973312	0.0163905	559	3.941118	4.005507	a
T1	12	3.980963	0.0163905	559	3.948768	4.013157	a
T3	12	4.004698	0.0163905	559	3.972503	4.036892	a
T2	15	3.994317	0.0163905	559	3.962123	4.026512	a
T1	15	4.003509	0.0163905	559	3.971315	4.035704	a
T3	15	4.030246	0.0163905	559	3.998051	4.062440	a
T2	18	4.015537	0.0163905	559	3.983343	4.047731	a
T1	18	4.028281	0.0163905	559	3.996086	4.060475	a
T3	18	4.059404	0.0163905	559	4.027209	4.091598	a
T2	21	4.079512	0.0163905	559	4.047317	4.111706	a
T3	21	4.097663	0.0163905	559	4.065469	4.129858	a
T1	21	4.097701	0.0163905	559	4.065506	4.129895	a
T3	26	4.191027	0.0163905	559	4.158833	4.223221	a
T2	26	4.197972	0.0163905	559	4.165777	4.230166	a
T1	26	4.230361	0.0163905	559	4.198167	4.262556	a
T2	33	4.279586	0.0163905	559	4.247392	4.311781	a
T1	33	4.321250	0.0163905	559	4.289056	4.353445	a
T3	33	4.326201	0.0163905	559	4.294007	4.358396	a
T2	36	4.364414	0.0163905	559	4.332219	4.396608	a
T3	36	4.390723	0.0163905	559	4.358528	4.422917	a
T1	36	4.404045	0.0163905	559	4.371851	4.436240	a
T3	40	4.407026	0.0163905	559	4.374832	4.439221	a
T2	40	4.432813	0.0163905	559	4.400619	4.465007	ab
T1	40	4.474017	0.0163905	559	4.441823	4.506212	b

within_means_tukey_Plot <- ggplot(dat, aes(y = weight, x = trt)) +
  facet_wrap(~ day, labeller = label_both) + # facet per G level
  geom_point(data = dat, aes(y = weight, x = trt, color = trt),
             alpha = 0.5, position = position_jitter(width = 0.1)) + # dots representing the raw data
  geom_point(data = within_means_tukey,aes(y = exp(weight), x = trt),
             color = "red", position = position_nudge(x = 0.2)) + # red dots representing the adjusted means
  geom_errorbar(data = within_means_tukey, aes(ymin = exp(lower.CL), ymax = exp(upper.CL), x = trt), color = "red", width = 0.1, position = position_nudge(x = 0.2)) + # red error bars representing the confidence limits of the adjusted means
  geom_text(data = within_means_tukey, aes(y = exp(weight), x = trt, label = str_trim(.group)), color = "black", position = position_nudge(x = 0.3), hjust = 0) + # red letters 
  scale_y_continuous(name = "Weight (g)", limits = c(30, NA), expand = expansion(mult = c(0, 0.1))) +
  scale_x_discrete(name = NULL) +
  theme_classic() + # clearer plot format 
  labs(caption = str_wrap("Opaque dots represent raw data. Red dots and error bars represent estimated marginal means ± 95% confidence interval per group. Means not sharing any letter are significantly different by the t-test at the 5% level of significance.", width = 120)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), legend.position = "bottom")
within_means_tukey_Plot

Experimental Design II

Conor Fair

2025-10-6

1 Experimental Design cont.

1.1 Split-plot Design

1.1.1 Review the Experimental Design

1.1.2 Importing the Data

1.1.3 Data Description

1.1.4 Exploring the Data

1.1.5 Fitting ANOVA Models with R

1.1.6 Sums of Squares Review

1.1.7 Assumptions

1.1.8 Mean Comparisons

1.1.9 Data Visualization

1.2 Repeated Measures

1.2.1 Review the Experimental Design

1.2.2 Importing the Data

1.2.3 Exploring the Data

1.2.4 Fitting ANOVA Models with R

1.2.5 Autocorrelated Errors

1.2.6 AR(1)

1.2.7 Compound Symmetry

1.2.8 Toeplitz

1.2.9 Unstructured Variance

1.2.10 Model Selection

1.2.11 Assumptions

1.2.12 Mean Comparisons

1.2.13 Data Visualization

1.3 Practice Example

1.3.1 Tips if you get stuck